// copyright:no copyright,use it if you like
// license boilerplate:BSD 3-Clause
// filename:Interpolation.hpp
// author:Tempest(https://gitee.com/Tempest56890)
// date:2022-10-20
// brief:implement the Newton formula and Hermite fomula.

#ifndef __TEMPEST_SPLINE__
#define __TEMPEST_SPLINE__

#include "Matrix.h"
#include " Interpolation.h"

class ppform : public Function {
private:
    std::vector<double> x,y;
    int N;
    std::vector<Polynomial> piecewise;
    double s_a,s_b,s_a_d,s_a_dd,s_b_d,s_b_dd;
public:
    ppform(std::vector<double> x,std::vector<double> y,double s_a=0.0,double s_b=0.0,double s_a_d=0.0,double s_a_dd=0.0,double s_b_d=0.0,double s_b_dd=0.0):x(x),y(y),s_a(s_a),s_b(s_b),s_a_d(s_a_d),s_a_dd(s_a_dd),s_b_d(s_b_d),s_b_dd(s_b_dd){
        if((int)x.size()<(int)y.size()){
            double ave=(x.back()-x.front())/(x.size());
            for(int i=1;i<=(int)y.size()-(int)x.size();i++){
                double x_end=x.back();
                x.push_back(x_end+ave);
            }
        }else if((int)x.size()>(int)y.size()){
            for(int i=1;i<=(int)x.size()-(int)y.size();i++){
                x.pop_back();
            }
        }
        s_a=y[0];s_b=y[y.size()-1];
        N=(int)x.size();
    };
    void print_linear_spline(){
        for(int i=0;i<=piecewise.size()-1;i++){
            std::cout<<"fl"<<(i+1)<<"="<<piecewise[i];
        }
    }
    void print_complete_spline(){
        for(int i=0;i<=piecewise.size()-1;i++){
            std::cout<<"fc"<<(i+1)<<"="<<piecewise[i];
        }
    }
    void print_natural_spline(){
        for(int i=0;i<=piecewise.size()-1;i++){
            std::cout<<"fn"<<(i+1)<<"="<<piecewise[i];
        }
    }
    void print_periodic_spline(){
        for(int i=0;i<=piecewise.size()-1;i++){
            std::cout<<"fp"<<(i+1)<<"="<<piecewise[i];
        }
    }
    void Linear_spline(){
        std::vector<Polynomial> temp;
        for(int i=0;i<=N-2;i++){
            temp.push_back(Polynomial(y[i]-x[i]*(y[i+1]-y[i])/(x[i+1]-x[i]),(y[i+1]-y[i])/(x[i+1]-x[i])));
        }
        piecewise=temp;
    }
    void Complete_spline(double ld=0.0,double rd=0.0,double ldd=0.0,double rdd=0.0){
        std::vector<Polynomial> temp;
        s_a_d=ld;s_b_d=rd;s_a_dd=ldd;s_b_dd=rdd;
        Matrix<double> A1(N-2,N-2),b1(N-2,1),sol1(N-2,1);
        A1.reset(0,0,2.0);A1.reset(0,1,(x[1]-x[0])/(x[2]-x[0]));
        A1.reset(N-3,N-4,(x[N-1]-x[N-2])/(x[N-1]-x[N-3]));A1.reset(N-3,N-3,2.0);
        b1.reset(0,0,3*((x[1]-x[0])/(x[2]-x[0]))*(y[2]-y[1])/(x[2]-x[1])+3*((x[2]-x[1])/(x[2]-x[0]))*(y[1]-y[0])/(x[1]-x[0])-s_a_d*(x[2]-x[1])/(x[2]-x[0]));
        b1.reset(N-3,0,3*((x[N-2]-x[N-3])/(x[N-1]-x[N-3]))*(y[N-1]-y[N-2])/(x[N-1]-x[N-2])+3*((x[N-1]-x[N-2])/(x[N-1]-x[N-3]))*(y[N-2]-y[N-3])/(x[N-2]-x[N-3])-s_b_d*(x[N-2]-x[N-3])/(x[N-1]-x[N-3]));
        for(int i=3;i<=N-2;i++){
            A1.reset(i-2,i-3,(x[i]-x[i-1])/(x[i]-x[i-2]));
            A1.reset(i-2,i-2,2.0);
            A1.reset(i-2,i-1,(x[i-1]-x[i-2])/(x[i]-x[i-2]));
            b1.reset(i-2,0,3*((x[i-2]-x[i-3])/(x[i-1]-x[i-3]))*(y[i-1]-y[i-2])/(x[i-1]-x[i-2])+3*((x[i-1]-x[i-2])/(x[i-1]-x[i-3]))*(y[i-2]-y[i-3])/(x[i-2]-x[i-3]));
        }
        sol1=A1.Matrix_inv(A1)*b1;
        Matrix<double> A2(N,N),b2(N,1),sol2(N,1);
        A2.reset(0,0,1.0);A2.reset(1,0,2.0);A2.reset(1,1,1.0);
        A2.reset(N-1,N-1,1.0);A2.reset(N-2,N-2,1.0);A2.reset(N-2,N-1,2.0);
        b2.reset(0,0,s_a_dd);b2.reset(N-1,0,s_b_dd);
        b2.reset(1,0,6*((y[1]-y[0])/(x[1]-x[0])-s_a_d)/(x[1]-x[0]));
        b2.reset(N-2,0,6*(s_b_d-(x[N-1]-x[N-2])/(x[N-1]-x[N-2]))/(x[N-1]-x[N-2]));
        for(int i=3;i<=N-2;i++){
            A2.reset(i-1,i,(x[i]-x[i-1])/(x[i]-x[i-2]));
            A2.reset(i-1,i-1,2.0);
            A2.reset(i-1,i-2,(x[i-1]-x[i-2])/(x[i]-x[i-2]));
            b2.reset(i-1,0,6*((y[i]-y[i-1])/(x[i]-x[i-1])-(y[i-1]-y[i-2])/(x[i-1]-x[i-2]))/(x[i]-x[i-2]));
        }
        sol2=A2.Matrix_inv(A2)*b2;
        std::vector<double> coef;
        coef.push_back(y[0]-s_a_d*x[0]+sol2(0,0)*x[0]*x[0]/2-pow(x[0],3)*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        coef.push_back(s_a_d-sol2(0,0)*x[0]+3*pow(x[0],2)*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        coef.push_back(sol2(0,0)/2-3*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3))*x[0]);
        coef.push_back((y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        temp.push_back(Polynomial(coef,3));
        coef.clear();
        for(int i=2;i<=N-1;i++){
            std::vector<double> coef;
            coef.push_back(y[i-1]-sol1(i-2,0)*x[i-1]+pow(x[i-1],2)*sol2(i-1,0)/2-pow(x[i-1],3)*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back(sol1(i-2,0)-sol2(i-1,0)*x[i-1]+3*pow(x[i-1],2)*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back(sol2(i-1,0)/2-3*x[i-1]*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back((y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            temp.push_back(Polynomial(coef,3));
            coef.clear();
        }
        piecewise=temp;
    }
    void Natural_spline(double ld=0.0,double rd=0.0,double ldd=0.0,double rdd=0.0){
        std::vector<Polynomial> temp;
        s_a_d=ld;s_b_d=rd;s_a_dd=0.0;s_b_dd=0.0;
        Matrix<double> A1(N-2,N-2),b1(N-2,1),sol1(N-2,1);
        A1.reset(0,0,2.0);A1.reset(0,1,(x[1]-x[0])/(x[2]-x[0]));
        A1.reset(N-3,N-4,(x[N-1]-x[N-2])/(x[N-1]-x[N-3]));A1.reset(N-3,N-3,2.0);
        b1.reset(0,0,3*((x[1]-x[0])/(x[2]-x[0]))*(y[2]-y[1])/(x[2]-x[1])+3*((x[2]-x[1])/(x[2]-x[0]))*(y[1]-y[0])/(x[1]-x[0])-s_a_d*(x[2]-x[1])/(x[2]-x[0]));
        b1.reset(N-3,0,3*((x[N-2]-x[N-3])/(x[N-1]-x[N-3]))*(y[N-1]-y[N-2])/(x[N-1]-x[N-2])+3*((x[N-1]-x[N-2])/(x[N-1]-x[N-3]))*(y[N-2]-y[N-3])/(x[N-2]-x[N-3])-s_b_d*(x[N-2]-x[N-3])/(x[N-1]-x[N-3]));
        for(int i=3;i<=N-2;i++){
            A1.reset(i-2,i-3,(x[i]-x[i-1])/(x[i]-x[i-2]));
            A1.reset(i-2,i-2,2.0);
            A1.reset(i-2,i-1,(x[i-1]-x[i-2])/(x[i]-x[i-2]));
            b1.reset(i-2,0,3*((x[i-2]-x[i-3])/(x[i-1]-x[i-3]))*(y[i-1]-y[i-2])/(x[i-1]-x[i-2])+3*((x[i-1]-x[i-2])/(x[i-1]-x[i-3]))*(y[i-2]-y[i-3])/(x[i-2]-x[i-3]));
        }
        sol1=A1.Matrix_inv(A1)*b1;
        Matrix<double> A2(N,N),b2(N,1),sol2(N,1);
        A2.reset(0,0,1.0);A2.reset(1,0,2.0);A2.reset(1,1,1.0);
        A2.reset(N-1,N-1,1.0);A2.reset(N-2,N-2,1.0);A2.reset(N-2,N-1,2.0);
        b2.reset(0,0,s_a_dd);b2.reset(N-1,0,s_b_dd);
        b2.reset(1,0,6*((y[1]-y[0])/(x[1]-x[0])-s_a_d)/(x[1]-x[0]));
        b2.reset(N-2,0,6*(s_b_d-(x[N-1]-x[N-2])/(x[N-1]-x[N-2]))/(x[N-1]-x[N-2]));
        for(int i=3;i<=N-2;i++){
            A2.reset(i-1,i,(x[i]-x[i-1])/(x[i]-x[i-2]));
            A2.reset(i-1,i-1,2.0);
            A2.reset(i-1,i-2,(x[i-1]-x[i-2])/(x[i]-x[i-2]));
            b2.reset(i-1,0,6*((y[i]-y[i-1])/(x[i]-x[i-1])-(y[i-1]-y[i-2])/(x[i-1]-x[i-2]))/(x[i]-x[i-2]));
        }
        sol2=A2.Matrix_inv(A2)*b2;
        std::vector<double> coef;
        coef.push_back(y[0]-s_a_d*x[0]+sol2(0,0)*x[0]*x[0]/2-pow(x[0],3)*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        coef.push_back(s_a_d-sol2(0,0)*x[0]+3*pow(x[0],2)*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        coef.push_back(sol2(0,0)/2-3*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3))*x[0]);
        coef.push_back((y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        temp.push_back(Polynomial(coef,3));
        coef.clear();
        for(int i=2;i<=N-1;i++){
            std::vector<double> coef;
            coef.push_back(y[i-1]-sol1(i-2,0)*x[i-1]+pow(x[i-1],2)*sol2(i-1,0)/2-pow(x[i-1],3)*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back(sol1(i-2,0)-sol2(i-1,0)*x[i-1]+3*pow(x[i-1],2)*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back(sol2(i-1,0)/2-3*x[i-1]*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back((y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            temp.push_back(Polynomial(coef,3));
            coef.clear();
        }
        piecewise=temp;
    }
    void Periodic_spline(double ld=0.0,double rd=0.0,double ldd=0.0,double rdd=0.0){
        std::vector<Polynomial> temp;
        s_b=s_a;
        s_a_d=ld;s_b_d=ld;s_a_dd=ldd;s_b_dd=ldd;
        Matrix<double> A1(N-2,N-2),b1(N-2,1),sol1(N-2,1);
        A1.reset(0,0,2.0);A1.reset(0,1,(x[1]-x[0])/(x[2]-x[0]));
        A1.reset(N-3,N-4,(x[N-1]-x[N-2])/(x[N-1]-x[N-3]));A1.reset(N-3,N-3,2.0);
        b1.reset(0,0,3*((x[1]-x[0])/(x[2]-x[0]))*(y[2]-y[1])/(x[2]-x[1])+3*((x[2]-x[1])/(x[2]-x[0]))*(y[1]-y[0])/(x[1]-x[0])-s_a_d*(x[2]-x[1])/(x[2]-x[0]));
        b1.reset(N-3,0,3*((x[N-2]-x[N-3])/(x[N-1]-x[N-3]))*(y[N-1]-y[N-2])/(x[N-1]-x[N-2])+3*((x[N-1]-x[N-2])/(x[N-1]-x[N-3]))*(y[N-2]-y[N-3])/(x[N-2]-x[N-3])-s_b_d*(x[N-2]-x[N-3])/(x[N-1]-x[N-3]));
        for(int i=3;i<=N-2;i++){
            A1.reset(i-2,i-3,(x[i]-x[i-1])/(x[i]-x[i-2]));
            A1.reset(i-2,i-2,2.0);
            A1.reset(i-2,i-1,(x[i-1]-x[i-2])/(x[i]-x[i-2]));
            b1.reset(i-2,0,3*((x[i-2]-x[i-3])/(x[i-1]-x[i-3]))*(y[i-1]-y[i-2])/(x[i-1]-x[i-2])+3*((x[i-1]-x[i-2])/(x[i-1]-x[i-3]))*(y[i-2]-y[i-3])/(x[i-2]-x[i-3]));
        }
        sol1=A1.Matrix_inv(A1)*b1;
        Matrix<double> A2(N,N),b2(N,1),sol2(N,1);
        A2.reset(0,0,1.0);A2.reset(1,0,2.0);A2.reset(1,1,1.0);
        A2.reset(N-1,N-1,1.0);A2.reset(N-2,N-2,1.0);A2.reset(N-2,N-1,2.0);
        b2.reset(0,0,s_a_dd);b2.reset(N-1,0,s_b_dd);
        b2.reset(1,0,6*((y[1]-y[0])/(x[1]-x[0])-s_a_d)/(x[1]-x[0]));
        b2.reset(N-2,0,6*(s_b_d-(x[N-1]-x[N-2])/(x[N-1]-x[N-2]))/(x[N-1]-x[N-2]));
        for(int i=3;i<=N-2;i++){
            A2.reset(i-1,i,(x[i]-x[i-1])/(x[i]-x[i-2]));
            A2.reset(i-1,i-1,2.0);
            A2.reset(i-1,i-2,(x[i-1]-x[i-2])/(x[i]-x[i-2]));
            b2.reset(i-1,0,6*((y[i]-y[i-1])/(x[i]-x[i-1])-(y[i-1]-y[i-2])/(x[i-1]-x[i-2]))/(x[i]-x[i-2]));
        }
        sol2=A2.Matrix_inv(A2)*b2;
        std::vector<double> coef;
        coef.push_back(y[0]-s_a_d*x[0]+sol2(0,0)*x[0]*x[0]/2-pow(x[0],3)*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        coef.push_back(s_a_d-sol2(0,0)*x[0]+3*pow(x[0],2)*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        coef.push_back(sol2(0,0)/2-3*(y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3))*x[0]);
        coef.push_back((y[1]-y[0]-s_a_d*(x[1]-x[0])-sol2(0,0)*pow(x[1]-x[0],2)/2)/(pow(x[1]-x[0],3)));
        temp.push_back(Polynomial(coef,3));
        coef.clear();
        for(int i=2;i<=N-1;i++){
            std::vector<double> coef;
            coef.push_back(y[i-1]-sol1(i-2,0)*x[i-1]+pow(x[i-1],2)*sol2(i-1,0)/2-pow(x[i-1],3)*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back(sol1(i-2,0)-sol2(i-1,0)*x[i-1]+3*pow(x[i-1],2)*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back(sol2(i-1,0)/2-3*x[i-1]*(y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            coef.push_back((y[i]-y[i-1]-sol1(i-2,0)*(x[i]-x[i-1])-sol2(i-1,0)*pow(x[i]-x[i-1],2)/2)/(pow(x[i]-x[i-1],3)));
            temp.push_back(Polynomial(coef,3));
            coef.clear();
        }
        piecewise=temp;
    }
    double operator()(const double t) const{
        if(t<x[0]||t>x[N-1]){
            std::cerr<<"Illegal input."<<std::endl;
            exit(-1);
        }else {
            int pos=0;
            while(t>x[pos+1]) pos++;
            return piecewise[pos](t);
        }
    }
};

class Bspline : public Function {
private:
    std::vector<double> x,y;
    int N;
    double x_1,x0,xn1,xn2,xn3;
    std::vector<Polynomial> B1,B2,B3;
    std::vector<Polynomial> piecewise;
    double s_a,s_b,s_a_d,s_a_dd,s_b_d,s_b_dd;
    std::vector<double> ai_coef;
public:
    Bspline(std::vector<double> x,std::vector<double> y,double s_a=0.0,double s_b=0.0,double s_a_d=0.0,double s_a_dd=0.0,double s_b_d=0.0,double s_b_dd=0.0):x(x),y(y),s_a(s_a),s_b(s_b),s_a_d(s_a_d),s_a_dd(s_a_dd),s_b_d(s_b_d),s_b_dd(s_b_dd){
        double ave=(x.back()-x.front())/(x.size());
        if((int)x.size()<(int)y.size()){
            for(int i=1;i<=(int)y.size()-(int)x.size();i++){
                double x_end=x.back();
                x.push_back(x_end+ave);
            }
        }else if((int)x.size()>(int)y.size()){
            for(int i=1;i<=(int)x.size()-(int)y.size();i++){
                x.pop_back();
            }
        }
        s_a=y[0];s_b=y[y.size()-1];
        N=(int)x.size();
        x0=x.front()-ave;x_1=x0-ave;
        xn1=x.back()+ave;xn2=xn1+ave;xn3=xn2+ave;
        B_base();
    };
    void B_base(){
        B1.push_back(Polynomial(-x_1/(x0-x_1),1/(x0-x_1)));
        B1.push_back(Polynomial(x[0]/(x[0]-x0),-1/(x[0]-x0)));
        B1.push_back(Polynomial(-x0/(x[0]-x0),1/(x[0]-x0)));
        B1.push_back(Polynomial(x[1]/(x[1]-x[0]),-1/(x[1]-x[0])));
        for(int i=2;i<=N-1;i++){
            B1.push_back(Polynomial(-x[i-2]/(x[i-1]-x[i-2]),1/(x[i-1]-x[i-2])));
            B1.push_back(Polynomial(x[i]/(x[i]-x[i-1]),-1/(x[i]-x[i-1])));
        }
        B1.push_back(Polynomial(-x[N-2]/(x[N-1]-x[N-2]),1/(x[N-1]-x[N-2])));
        B1.push_back(Polynomial(xn1/(xn1-x[N-1]),-1/(xn1-x[N-1])));
        B1.push_back(Polynomial(-x[N-1]/(xn1-x[N-1]),1/(xn1-x[N-1])));
        B1.push_back(Polynomial(xn2/(xn2-xn1),-1/(xn2-xn1)));
        B1.push_back(Polynomial(-xn1/(xn2-xn1),1/(xn2-xn1)));
        B1.push_back(Polynomial(xn3/(xn3-xn2),-1/(xn3-xn2)));
        //
        std::vector<double> temp;
        temp.push_back(x_1);temp.push_back(x0);
        for(int i=0;i<=N-1;i++) temp.push_back(x[i]);
        temp.push_back(xn1);temp.push_back(xn2);temp.push_back(xn3);
        for(int i=1;i<=N+1;i++){
            B2.push_back(Polynomial(pow(temp[i-1],2)/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1])),-2*temp[i-1]/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1])),1/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1]))));
            B2.push_back(Polynomial(-temp[i+1]*temp[i-1]/((temp[i+1]-temp[i-1])*(temp[i+1]-temp[i]))-temp[i]*temp[i+2]/((temp[i+2]-temp[i])*(temp[i+1]-temp[i])),2/(temp[i+1]-temp[i]),-1/((temp[i+1]-temp[i-1])*(temp[i+1]-temp[i]))-1/((temp[i+2]-temp[i])*(temp[i+1]-temp[i]))));
            B2.push_back(Polynomial(pow(temp[i+2],2)/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1])),-2*temp[i+2]/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1])),1/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1]))));
        }
    }
    double b_operator(const double t) {
        std::vector<double> temp;
        temp.push_back(x_1);temp.push_back(x0);
        for(int i=0;i<=N-1;i++) temp.push_back(x[i]);
        temp.push_back(xn1);temp.push_back(xn2);temp.push_back(xn3);
        if(t<temp[0]||t>temp[N+3]){
            std::cerr<<"Illegal input."<<std::endl;
            exit(-1);
        }else if(t<temp[1]){
            return (B2[0](t));
        }else if(t<temp[2]){
            return (B2[1](t)+B2[3](t));
        }else if(t>temp[N+2]){
            return (B2[3*N+2](t));
        }else if(t>temp[N+1]){
            return (B2[3*N+1](t)+B2[3*N-1](t));
        }else {
            int pos=2;
            while(t>=temp[pos+1]) pos++;
            pos--;
            return (B2[3*pos-1](t)+B2[3*pos+1](t)+B2[3*pos+3](t));
        }
    }
    void Complete_spline(double ld=0.0,double rd=0.0,double ldd=0.0,double rdd=0.0){
        s_a_d=ld;
        Matrix<double> A(N+1,N+1),b(N+1,1),sol(N+1,1);
        A.reset(1,0,B2[1](x[0]));A.reset(1,1,B2[3](x[0]));b.reset(0,0,s_a_d);
        A.reset(0,0,-(x[0]+x_1)/((x[0]-x_1)*(x[0]-x0))+(x[1]-2*x[0]-x0)/((x[1]-x0)*(x[0]-x0)));
        A.reset(0,1,2/(x[1]-x0));
        for(int i=2;i<=N;i++){
            A.reset(i,i-2,B2[3*i-4](x[i-1]));
            A.reset(i,i-1,B2[3*i-2](x[i-1]));
            A.reset(i,i,B2[3*i](x[i-1]));
        }
        for(int i=1;i<=N;i++) b.reset(i,0,y[i-1]);
        sol=A.Matrix_inv(A)*b;
        for(int i=0;i<=N;i++){
            ai_coef.push_back(sol(i,0));
        }
        for(int i=1;i<=N+3;i++) {
            piecewise.push_back(Polynomial(0.0));
        }
        for(int i=0;i<=N;i++){
            piecewise[i]=piecewise[i]+B2[3*i];
            piecewise[i+1]=piecewise[i]+B2[3*i+1];
            piecewise[i+2]=piecewise[i]+B2[3*i+2];
        }
    }
    void Natural_spline(double ld=0.0,double rd=0.0,double ldd=0.0,double rdd=0.0){
        s_a_d=ld;s_b_d=rd;s_a_dd=0.0;s_b_dd=0.0;
        Matrix<double> A(N+1,N+1),b(N+1,1),sol(N+1,1);
        A.reset(1,0,B2[1](x[0]));A.reset(1,1,B2[3](x[0]));b.reset(0,0,s_a_d);
        A.reset(0,0,-(x[0]+x_1)/((x[0]-x_1)*(x[0]-x0))+(x[1]-2*x[0]-x0)/((x[1]-x0)*(x[0]-x0)));
        A.reset(0,1,2/(x[1]-x0));
        for(int i=2;i<=N;i++){
            A.reset(i,i-2,B2[3*i-4](x[i-1]));
            A.reset(i,i-1,B2[3*i-2](x[i-1]));
            A.reset(i,i,B2[3*i](x[i-1]));
        }
        for(int i=1;i<=N;i++) b.reset(i,0,y[i-1]);
        sol=A.Matrix_inv(A)*b;
        for(int i=0;i<=N;i++){
            ai_coef.push_back(sol(i,0));
        }
        for(int i=1;i<=N+3;i++) {
            piecewise.push_back(Polynomial(0.0));
        }
        for(int i=0;i<=N;i++){
            piecewise[i]=piecewise[i]+B2[3*i];
            piecewise[i+1]=piecewise[i]+B2[3*i+1];
            piecewise[i+2]=piecewise[i]+B2[3*i+2];
        }
    }
    void Periodic_spline(double ld=0.0,double rd=0.0,double ldd=0.0,double rdd=0.0){
        s_b=s_a;
        s_a_d=ld;s_b_d=ld;s_a_dd=ldd;s_b_dd=ldd;
        Matrix<double> A(N+1,N+1),b(N+1,1),sol(N+1,1);
        A.reset(1,0,B2[1](x[0]));A.reset(1,1,B2[3](x[0]));b.reset(0,0,s_a_d);
        A.reset(0,0,-(x[0]+x_1)/((x[0]-x_1)*(x[0]-x0))+(x[1]-2*x[0]-x0)/((x[1]-x0)*(x[0]-x0)));
        A.reset(0,1,2/(x[1]-x0));
        for(int i=2;i<=N;i++){
            A.reset(i,i-2,B2[3*i-4](x[i-1]));
            A.reset(i,i-1,B2[3*i-2](x[i-1]));
            A.reset(i,i,B2[3*i](x[i-1]));
        }
        for(int i=1;i<=N;i++) b.reset(i,0,y[i-1]);
        sol=A.Matrix_inv(A)*b;
        for(int i=0;i<=N;i++){
            ai_coef.push_back(sol(i,0));
        }
        for(int i=1;i<=N+3;i++) {
            piecewise.push_back(Polynomial(0.0));
        }
        for(int i=0;i<=N;i++){
            piecewise[i]=piecewise[i]+B2[3*i];
            piecewise[i+1]=piecewise[i]+B2[3*i+1];
            piecewise[i+2]=piecewise[i]+B2[3*i+2];
        }
    }
    double operator()(const double t) const{
        if(t<x_1||t>xn2){
            std::cerr<<"Illegal input."<<std::endl;
            exit(-1);
        }else {
            std::vector<double> temp;
            temp.push_back(x_1);temp.push_back(x0);
            for(int i=1;i<=N;i++) temp.push_back(x[i-1]);
            ppform piecewice(x,y);piecewice.Complete_spline();
            temp.push_back(xn1);temp.push_back(xn2);temp.push_back(xn3);
            int pos=0;
            while(t>temp[pos+1]) pos++;
            return piecewice(t);
        }
    }
};

class Curve:public Function
{
private:
    std::vector<double> x1,y1,x2,y2;
    int N;
    double x_1,x0,xn1,xn2,xn3;
    std::vector<Polynomial> B1,B2,B3;
    std::vector<Polynomial> piecewise1,piecewise2;
    double s_a,s_b,s_a_d,s_a_dd,s_b_d,s_b_dd;
    std::vector<double> ai_coef;
    std::vector<double> characteristic_points;
public:
    Curve(std::vector<double> x1,std::vector<double> y1,std::vector<double> x2,std::vector<double> y2,std::vector<double> c_p,double s_a=0.0,double s_b=0.0,double s_a_d=0.0,double s_a_dd=0.0,double s_b_d=0.0,double s_b_dd=0.0):x1(x1),y1(y1),x2(x2),y2(y2),characteristic_points(c_p),s_a(s_a),s_b(s_b),s_a_d(s_a_d),s_a_dd(s_a_dd),s_b_d(s_b_d),s_b_dd(s_b_dd){
        double ave=(x1.back()-x1.front())/(x1.size());
        if((int)x1.size()<(int)y1.size()){
            for(int i=1;i<=(int)y1.size()-(int)x1.size();i++){
                double x1_end=x1.back();
                x1.push_back(x1_end+ave);
            }
        }else if((int)x1.size()>(int)y1.size()){
            for(int i=1;i<=(int)x1.size()-(int)y1.size();i++){
                x1.pop_back();
            }
        }
        s_a=y1[0];s_b=y1[y1.size()-1];
        N=(int)x1.size();
        x0=x1.front()-ave;x_1=x0-ave;
        xn1=x1.back()+ave;xn2=xn1+ave;xn3=xn2+ave;
        Curve_base1();
        curve_base2();
    };
    void Curve_base1(){
        B1.push_back(Polynomial(-x_1/(x0-x_1),1/(x0-x_1)));
        B1.push_back(Polynomial(x1[0]/(x1[0]-x0),-1/(x1[0]-x0)));
        B1.push_back(Polynomial(-x0/(x1[0]-x0),1/(x1[0]-x0)));
        B1.push_back(Polynomial(x1[1]/(x1[1]-x1[0]),-1/(x1[1]-x1[0])));
        for(int i=2;i<=N-1;i++){
            B1.push_back(Polynomial(-x1[i-2]/(x1[i-1]-x1[i-2]),1/(x1[i-1]-x1[i-2])));
            B1.push_back(Polynomial(x1[i]/(x1[i]-x1[i-1]),-1/(x1[i]-x1[i-1])));
        }
        B1.push_back(Polynomial(-x1[N-2]/(x1[N-1]-x1[N-2]),1/(x1[N-1]-x1[N-2])));
        B1.push_back(Polynomial(xn1/(xn1-x1[N-1]),-1/(xn1-x1[N-1])));
        B1.push_back(Polynomial(-x1[N-1]/(xn1-x1[N-1]),1/(xn1-x1[N-1])));
        B1.push_back(Polynomial(xn2/(xn2-xn1),-1/(xn2-xn1)));
        B1.push_back(Polynomial(-xn1/(xn2-xn1),1/(xn2-xn1)));
        B1.push_back(Polynomial(xn3/(xn3-xn2),-1/(xn3-xn2)));
        //
        std::vector<double> temp;
        temp.push_back(x_1);temp.push_back(x0);
        for(int i=0;i<=N-1;i++) temp.push_back(x1[i]);
        temp.push_back(xn1);temp.push_back(xn2);temp.push_back(xn3);
        for(int i=1;i<=N+1;i++){
            B2.push_back(Polynomial(pow(temp[i-1],2)/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1])),-2*temp[i-1]/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1])),1/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1]))));
            B2.push_back(Polynomial(-temp[i+1]*temp[i-1]/((temp[i+1]-temp[i-1])*(temp[i+1]-temp[i]))-temp[i]*temp[i+2]/((temp[i+2]-temp[i])*(temp[i+1]-temp[i])),2/(temp[i+1]-temp[i]),-1/((temp[i+1]-temp[i-1])*(temp[i+1]-temp[i]))-1/((temp[i+2]-temp[i])*(temp[i+1]-temp[i]))));
            B2.push_back(Polynomial(pow(temp[i+2],2)/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1])),-2*temp[i+2]/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1])),1/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1]))));
        }
    }
    void curve_base2(){
        B1.push_back(Polynomial(-x_1/(x0-x_1),1/(x0-x_1)));
        B1.push_back(Polynomial(x2[0]/(x2[0]-x0),-1/(x2[0]-x0)));
        B1.push_back(Polynomial(-x0/(x2[0]-x0),1/(x2[0]-x0)));
        B1.push_back(Polynomial(x2[1]/(x2[1]-x2[0]),-1/(x2[1]-x2[0])));
        for(int i=2;i<=N-1;i++){
            B1.push_back(Polynomial(-x2[i-2]/(x2[i-1]-x2[i-2]),1/(x2[i-1]-x2[i-2])));
            B1.push_back(Polynomial(x2[i]/(x2[i]-x2[i-1]),-1/(x2[i]-x2[i-1])));
        }
        B1.push_back(Polynomial(-x2[N-2]/(x2[N-1]-x2[N-2]),1/(x2[N-1]-x2[N-2])));
        B1.push_back(Polynomial(xn1/(xn1-x2[N-1]),-1/(xn1-x2[N-1])));
        B1.push_back(Polynomial(-x2[N-1]/(xn1-x2[N-1]),1/(xn1-x2[N-1])));
        B1.push_back(Polynomial(xn2/(xn2-xn1),-1/(xn2-xn1)));
        B1.push_back(Polynomial(-xn1/(xn2-xn1),1/(xn2-xn1)));
        B1.push_back(Polynomial(xn3/(xn3-xn2),-1/(xn3-xn2)));
        //
        std::vector<double> temp;
        temp.push_back(x_1);temp.push_back(x0);
        for(int i=0;i<=N-1;i++) temp.push_back(x2[i]);
        temp.push_back(xn1);temp.push_back(xn2);temp.push_back(xn3);
        for(int i=1;i<=N+1;i++){
            B2.push_back(Polynomial(pow(temp[i-1],2)/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1])),-2*temp[i-1]/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1])),1/((temp[i+1]-temp[i-1])*(temp[i]-temp[i-1]))));
            B2.push_back(Polynomial(-temp[i+1]*temp[i-1]/((temp[i+1]-temp[i-1])*(temp[i+1]-temp[i]))-temp[i]*temp[i+2]/((temp[i+2]-temp[i])*(temp[i+1]-temp[i])),2/(temp[i+1]-temp[i]),-1/((temp[i+1]-temp[i-1])*(temp[i+1]-temp[i]))-1/((temp[i+2]-temp[i])*(temp[i+1]-temp[i]))));
            B2.push_back(Polynomial(pow(temp[i+2],2)/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1])),-2*temp[i+2]/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1])),1/((temp[i+2]-temp[i])*(temp[i+2]-temp[i+1]))));
        }
    }
    void Curve_spline(){
        std::vector<Polynomial> temp1,temp2;
        Matrix<double> A1(N-2,N-2),b1(N-2,1),sol1(N-2,1);
        A1.reset(0,0,2.0);A1.reset(0,1,(x1[1]-x1[0])/(x1[2]-x1[0]));
        A1.reset(N-3,N-4,(x1[N-1]-x1[N-2])/(x1[N-1]-x1[N-3]));A1.reset(N-3,N-3,2.0);
        b1.reset(0,0,3*((x1[1]-x1[0])/(x1[2]-x1[0]))*(y1[2]-y1[1])/(x1[2]-x1[1])+3*((x1[2]-x1[1])/(x1[2]-x1[0]))*(y1[1]-y1[0])/(x1[1]-x1[0])-s_a_d*(x1[2]-x1[1])/(x1[2]-x1[0]));
        b1.reset(N-3,0,3*((x1[N-2]-x1[N-3])/(x1[N-1]-x1[N-3]))*(y1[N-1]-y1[N-2])/(x1[N-1]-x1[N-2])+3*((x1[N-1]-x1[N-2])/(x1[N-1]-x1[N-3]))*(y1[N-2]-y1[N-3])/(x1[N-2]-x1[N-3])-s_b_d*(x1[N-2]-x1[N-3])/(x1[N-1]-x1[N-3]));
        for(int i=3;i<=N-2;i++){
            A1.reset(i-2,i-3,(x1[i]-x1[i-1])/(x1[i]-x1[i-2]));
            A1.reset(i-2,i-2,2.0);
            A1.reset(i-2,i-1,(x1[i-1]-x1[i-2])/(x1[i]-x1[i-2]));
            b1.reset(i-2,0,3*((x1[i-2]-x1[i-3])/(x1[i-1]-x1[i-3]))*(y1[i-1]-y1[i-2])/(x1[i-1]-x1[i-2])+3*((x1[i-1]-x1[i-2])/(x1[i-1]-x1[i-3]))*(y1[i-2]-y1[i-3])/(x1[i-2]-x1[i-3]));
        }
        sol1=A1.Matrix_inv(A1)*b1;
        Matrix<double> A2(N,N),b2(N,1),sol2(N,1);
        A2.reset(0,0,1.0);A2.reset(1,0,2.0);A2.reset(1,1,1.0);
        A2.reset(N-1,N-1,1.0);A2.reset(N-2,N-2,1.0);A2.reset(N-2,N-1,2.0);
        b2.reset(0,0,s_a_dd);b2.reset(N-1,0,s_b_dd);
        b2.reset(1,0,6*((y1[1]-y1[0])/(x1[1]-x1[0])-s_a_d)/(x1[1]-x1[0]));
        b2.reset(N-2,0,6*(s_b_d-(x1[N-1]-x1[N-2])/(x1[N-1]-x1[N-2]))/(x1[N-1]-x1[N-2]));
        for(int i=3;i<=N-2;i++){
            A2.reset(i-1,i,(x1[i]-x1[i-1])/(x1[i]-x1[i-2]));
            A2.reset(i-1,i-1,2.0);
            A2.reset(i-1,i-2,(x1[i-1]-x1[i-2])/(x1[i]-x1[i-2]));
            b2.reset(i-1,0,6*((y1[i]-y1[i-1])/(x1[i]-x1[i-1])-(y1[i-1]-y1[i-2])/(x1[i-1]-x1[i-2]))/(x1[i]-x1[i-2]));
        }
        sol2=A2.Matrix_inv(A2)*b2;
        std::vector<double> coef;
        coef.push_back(y2[0]-s_a_d*x2[0]+sol2(0,0)*x2[0]*x2[0]/2-pow(x2[0],3)*(y2[1]-y2[0]-s_a_d*(x2[1]-x2[0])-sol2(0,0)*pow(x2[1]-x2[0],2)/2)/(pow(x2[1]-x2[0],3)));
        coef.push_back(s_a_d-sol2(0,0)*x2[0]+3*pow(x2[0],2)*(y2[1]-y2[0]-s_a_d*(x2[1]-x2[0])-sol2(0,0)*pow(x2[1]-x2[0],2)/2)/(pow(x2[1]-x2[0],3)));
        coef.push_back(sol2(0,0)/2-3*(y2[1]-y2[0]-s_a_d*(x2[1]-x2[0])-sol2(0,0)*pow(x2[1]-x2[0],2)/2)/(pow(x2[1]-x2[0],3))*x2[0]);
        coef.push_back((y2[1]-y2[0]-s_a_d*(x2[1]-x2[0])-sol2(0,0)*pow(x2[1]-x2[0],2)/2)/(pow(x2[1]-x2[0],3)));
        temp1.push_back(Polynomial(coef,3));
        coef.clear();
        for(int i=2;i<=N-1;i++){
            std::vector<double> coef;
            coef.push_back(y2[i-1]-sol1(i-2,0)*x2[i-1]+pow(x2[i-1],2)*sol2(i-1,0)/2-pow(x2[i-1],3)*(y2[i]-y2[i-1]-sol1(i-2,0)*(x2[i]-x2[i-1])-sol2(i-1,0)*pow(x2[i]-x2[i-1],2)/2)/(pow(x2[i]-x2[i-1],3)));
            coef.push_back(sol1(i-2,0)-sol2(i-1,0)*x2[i-1]+3*pow(x2[i-1],2)*(y2[i]-y2[i-1]-sol1(i-2,0)*(x2[i]-x2[i-1])-sol2(i-1,0)*pow(x2[i]-x2[i-1],2)/2)/(pow(x2[i]-x2[i-1],3)));
            coef.push_back(sol2(i-1,0)/2-3*x2[i-1]*(y2[i]-y2[i-1]-sol1(i-2,0)*(x2[i]-x2[i-1])-sol2(i-1,0)*pow(x2[i]-x2[i-1],2)/2)/(pow(x2[i]-x2[i-1],3)));
            coef.push_back((y2[i]-y2[i-1]-sol1(i-2,0)*(x2[i]-x2[i-1])-sol2(i-1,0)*pow(x2[i]-x2[i-1],2)/2)/(pow(x2[i]-x2[i-1],3)));
            temp2.push_back(Polynomial(coef,3));
            coef.clear();
        }
        piecewise1=temp1;
        piecewise2=temp2;
    }
    double Curve_operator(const double t) {
        std::vector<double> temp;
        temp.push_back(x_1);temp.push_back(x0);
        for(int i=0;i<=N-1;i++) temp.push_back(x1[i]);
        temp.push_back(xn1);temp.push_back(xn2);temp.push_back(xn3);
        if(t<temp[0]||t>temp[N+3]){
            std::cerr<<"Illegal input."<<std::endl;
            exit(-1);
        }else if(t<temp[1]){
            return (B2[0](t));
        }else if(t<temp[2]){
            return (B2[1](t)+B2[3](t));
        }else if(t>temp[N+2]){
            return (B2[3*N+2](t));
        }else if(t>temp[N+1]){
            return (B2[3*N+1](t)+B2[3*N-1](t));
        }else {
            int pos=2;
            while(t>=temp[pos+1]) pos++;
            pos--;
            return (B2[3*pos-1](t)+B2[3*pos+1](t)+B2[3*pos+3](t));
        }
    }
};

#else
#endif /* Spline_h */
